A simulation study of localization of electromagnetic waves in two-dimensional 

random dipolar systems 
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We study the propagation and scattering of electromagnetic waves by random arrays of dipolar 
cylinders in a uniform medium. A set of self-consistent equations, incorporating all orders of multiple 
scattering of the electromagnetic waves, is derived from first principles and then solved numerically 
for electromagnetic fields. For certain ranges of frequencies, spatially localized electromagnetic 
waves appear in such a simple but realistic disordered system. Dependence of localization on the 
frequency, radiation damping, and filling factor is shown. The spatial behavior of the total, coherent 
and diffusive waves is explored in detailed, and found to be comply with a physical intuitive picture. 
A phase diagram characterizing localization is presented, in agreement with previous investigations 
on other systems. 

PACS numbers: 42.25.Hz, 41.90.+e 



I. INTRODUCTION 



The concept of localization was originally introduced 
by Anderson 1] for electrons in a crystal. In the case of a 
perfectly periodic lattice, except in the gaps all the elec- 
tronic states are extended and are represented by Bloch 
states. When a sufficient amount of disorders is added 
to the lattice, for example in the form of random poten- 
tials, the electrons may become spatially localized due to 
the multiple scattering by the disorders. In such a case, 
the eigenstates are exponentially confined in the spacejjj. 
The inception of the localization concept has opened a 
new era for the study of electrons in disordered systems, 
and stimulated a tremendous research. The concept of lo- 
calization has also rendered a great development in many 
other fields such as seismology Q, oceanologyQ, and ran- 
dom lasers^, to name just a few. The great efforts have 
been summarized in a number of excellent reviews (e. g. 

Over the past two decades, localization of classical 
waves has been under intensive investigations, lea ding 

.^BSj,IIllil 

miM Ha). Such a 
localization phenomenon has been characterized by two 
levels. One is the weak localization associated with the 
enhanced backscattering. That is, waves which propa- 
gate in the two opposite directions along a loop will ob- 
tain the same phase and interfere constructively at the 
emission site, thus enhancing the backscattering. The 
second is the strong localization, without confusion often 
just termed as localization, in which a significant inhi- 
bition of transmission appears and the energy is mostly 
confined spatially in the vicinity of the emission site. 
While the weak localization, regarded as a precursor to 



'Previously publishing as Kang-Xin Wang 
^Electronic address: zhen@phy. ncu.edu. tw| 



the strong localization, has been well studied both the- 
oretically (e. g. the monograph^) and experimentally 
(e. g. 25]), observation of strong localization of classical 
waves for hig her than one dimension remains a subject 
of debate jUES H3- primarily because a suitable sys- 
tem is hard to find and the observation is often obscured 
by such effects as the residual absorption (2^ and scat- 
tering attenuation. Asides from few exceptions [l^, most 
experiments are based on observations of the exponen- 
tial decay of waves as they propagate through disordered 
media. This was pointed out in Ref. j^j. According 
to |21j, this type of experiments is not sufficient to dis- 
cern whether the medium really only has localized states. 
Unwanted effects of non-localization origin can also con- 
tribute to the exponential decay, making data interpre- 
tation possibly ambiguous. In a conclusion, Sigalas et al. 
|22| pointed out that there is no conclusive experimental 
evidence for localization of electromagnetic waves (EM) 
in two dimensions (2D). We mention that there was a re- 
port of the observation of microwave localization in two 
dimensions when a transmitting source is inside disor- 
dered media |Tgj. 

The difficulties in observing localization of EM waves 
mainly lie in a couple of problems. First, wave local- 
ization only appears for strongly scattering media, and 
such a medium is often hard to find. Second, localiza- 
tion effects arc often entangled with other effects such 
as dissipation, wave deflection, or boundary effects[l2j|. 
making data interpretation often ambiguous. In a recent 
communication, a simple but seemingly realistic model 
system has been proposed to study EM localization in 2D 
random media[3(j. This model originated from the pre- 
vious study of the radiative effects of the electric dipoles 
embedded in structured cavities |3j. It was shown that 
EM localization is possible in such a disordered system. 
When localization occurs, a coherent behavior appears 
and is revealed as a unique property differentiating lo- 
calization from either the residual absorption or the at- 
tenuation effects. 
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With the present paper, we wish to explore further the 
localization properties of the system outlined in |3(j • We 
will investigate the dependence of localization behavior 
on a numbers of parameters including frequency, filling 
factor, scattering strength, damping effect, and two dif- 
ferent ways of measuring localization. Additionally, the 
spatial behaviors of the total, coherent and diffusive wave 
intensities will be studied, and are shown to comply with 
a simple physical intuition. 



the sake of convenience and completeness. 

Upon stimulation, each dipole will radiate EM waves. 
The radiated waves will then repeatedly interact with the 
dipoles, forming a process of multiple scattering. The 
equation of motion for the i-th dipole is 



dt 2 ' 



q 2 d 
—E z (fi) - b ,i—Pi 
m,- dt 



for i = 1,2,. 



(1) 



II. THE SYSTEM AND THEORETICAL 
FORMULATION 

A. The system 

Following Erdogan et al.j^j], we consider 2D dipoles 
as an ensemble of harmonically bound charge elements. 
In this way, each 2D dipole is regarded as a single dipole 
line, characterized by the charge and dipole moment per 
unit length. Assume that N parallel dipole lines, aligned 
along the z-axis, are embedded in a uniform dielectric 
medium and randomly located at = 1, 2, ... , N). The 
averaged distance between dipoles is d. A stimulating 
dipole line source is located at f s , transmitting a contin- 
uous wave of angular frequency to. By the geometrical 
symmetry of the system, we only need to consider the z 
component of the electrical waves. The conceptual layout 
of the system is illustrated by Fig. 

Here it is shown that a line source is located at the ori- 
gin, and the receiver is located at some distance from the 
source. After multiply scattered, the transmitted waves 
reach at the receiver. 
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where wo, t is the resonance (natural) frequency, pi, qi 
and rrii the dipole moment, charge and effective mass per 
unit length of the i-th dipole respectively. E z {r%) is the 
total electrical field acting on dipole pi, which includes 
the radiated field from other dipoles and also the directly 
field from the source. The factor &o,i denotes the damping 
due to energy loss and radiation, and can be determined 
by energy conservation. Without energy loss (to such as 
heat), &cm can be determined from the balance between 
the radiative and vibrational energies, and is given as|3lj 



bo,i 



(2) 



with e being the constant permittivity and c the speed of 
light in the medium separately. 

Equation is virtually the same as Eq. (1) in [3l| . 
The only difference is that in [3l| , E z is the reflected field 
at the dipole due to the presence of reflecting surround- 
ing structures, while in the present case the field is from 
the stimulating source and the radiation from all other 
dipoles. 

The transmitted electrical field from the continuous 
line source is determined by the Maxwell equations 



d 2 



c 2 dt 2 



j G (f-r s ) = -4,uow 2 po7rrJ (2 ) {r-r s )e 



(3) 

where u> is the radiation frequency, and po is the source 
strength and is set to be unit. The solution of Eq. is 
clearly 



G (f- f° s ) = (Li a uj 2 p a )iTtH { Q\k\r - f 



with k = ui/c, and Hq 1 ^ being the zero-th order Hankel 
function of the first kind. 

Similarly, the radiated field from the i-th dipole is 
given by 



(4) 



FIG. 1: Conceptual layout of the system and the simulation. 



v2 - ^) Gi{f- n) = ^ Pi 6W(r- n). (5) 



B. The formulation 



Although much of the following materials can be re- 
ferred to in [30j, we repeat the important parts here for 



The field arriving at the i-th dipole is composed of the 
direct field from the source and the radiation from all 
other dipoles, and thus is given as 



E z (fi) = G (n - r s ) + 



N 

E 



Gjin-rj). (6) 
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Substituting Eqs. (2J|, ©, an d © i nto Eq. ©, and 
writing = p^e~ , we arrive at 

(-w 2 + luI , - iujb ,i)Pi = 



N 



(7) 



This set of linear equations can be solved numerically for 
Pi. After pi are obtained, the total field at any space 
point can be readily calculated from 



N 



(8) 



3 = 1 



In the standard approach to wave localization, waves 
are said to be localized when the square modulus of the 
field \E(r)\ 2 , representing the wave energy, is spatially 
localized after the trivial cylindrically spreading effect is 
eliminated. Obviously, this is equivalent to say that the 
further away is the dipole from the source, the smaller its 
oscillation amplitude, expected to follow an exponentially 
decreasing pattern. 

To this end, it is instructive to point out that an alter- 
native two dimensional dipole model was devised previ- 
ously by Rusek and Orlowski 21j. The authors derived a 
set of linear algebraic equations, which is similar in form 
to the above Eq. ©. However, there are some funda- 
mental discrepancies between the two models. In |2lj . 
the interaction between dipoles and the external field is 
derived by the energy conservation, while in the present 
case the coupling is determined without ambiguity by the 
Newton's second law. The former leads to an undeter- 
mined phase factor. According to, e. g. Refs. [3ll 13^. 
the energy conservation can only give the radiation fac- 
tor in Eq. ©. We would also like to point out that 
the set of couple equations in Eq. © is similar in spirit 
to the tight-binding model used to study the electronic 
localization^ l33|. 

There are several ways to introduce randomness to 
Eq. ©. For example, the disorder may be introduced by 
randomly varying such properties of individual dipoles as 
the charge, the mass or the two combined. This is the 
most common way that the disorder is introduced into 
the tight-binding model for electronic waves |33j. In the 
present study, the disorder is brought in by the random 
distribution of the dipoles. 

For simplicity yet without losing generality, assume 
that all the dipoles are identical and they are randomly 
distributed within a square area. The source is located at 
the center (set to be the origin) of this area. For conve- 
nience, we make Eq. © non-dimensional by scaling the 
frequency by the resonance frequency of a single dipole 
luq. This will lead to two independent non-dimensional 

parameters b — and b' = ^ f^)' Both parame- 
ters may be adjusted in the experiment. For example, the 
factor bo can be modified by coating layered structures 



around the dipoles[3lJ. Then Eq. © becomes simply 

(-/ 2 + 1 - ib' ) Pi = 



ibf 



N 



PoH^ikln - r s \) + PiH^im-fj]) 



(9) 



with / = Eq. © is self-consistent and can be solved 
numerically for pi and then the total field is obtained 
through Eqs. ©, © and ©. 



III. THE RESULTS AND DISCUSSION 
A. Two numerical measuring scenarios 

In the following computation, we will consider two sce- 
narios. They are illustrated in Fig.[2with the coordinates 
being shown. In both cases, the sample takes a fixed rect- 
angular shape of which the size may vary. The dipoles, 
denoted by the small circles, are placed within the rect- 
angle in a complete randomness. The receiver, denoted 
by the filled black circle, is placed on the x-axis. 

In case (A), the sample size is fixed. The receiver is 
placed along x-axis to measure the signal at various posi- 
tions, yielding the result of the transmission signal versus 
the distance between the source and the receiver, termed 
as the travelling distance in the paper. This scenario 
complies with the original conjecture of observation of 
localization^. In case (B), the receiver is placed at a 
very small distance outside the sample. While the sam- 
ple size varies, this method is to measure the transmission 
through the sample of various sizes, yielding the result of 
the transmission signal versus the sample size. 
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FIG. 2: The conceptual layout of the two measuring methods. 



B. General discussion of localization 

Before moving to solve Eq. © for the phenomenon of 
localization of EM waves, we discuss some general prop- 
erties of wave localization. 

The coherence in localization Although some parts of 
the following discussion has been reported earlier, we re- 
peat here for the sake of convenience and importance. 



4 



The energy flow of EM waves is J ~ E x H . By invoking 
the Maxwell equations to relate the electrical and mag- 
netic fields, we can derive that the time averaged energy 
flow is 



< J >t 



dtJ ~ \E\ 2 ve, 



(10) 



where the electrical field is written as E = eE\E\e , with 
&e denoting the direction, \E\ and 9 being the amplitude 
and the phase respectively. It is clear from Eq. Q10[l that 
when 9 is constant, at least by spatial domains, while 
\E\ ^ 0, the flow would come to a stop and the energy 
will be localized or stored in the space. In the localized 
state, a source can no longer radiate energies. Alter- 
natively, we can write the oscillation of the dipoles as 
Pi = \Pi\e l6i - By studying the square modulus of Pi in 
the form of — r s ||pj| 2 , and its phase 6^, we can also 
investigate the localization of EM waves. Note here that 
the factor |r*j — r s \ is to eliminate the cylindrical spread- 
ing effect in 2D as can be seen from the expansion of 
the Hankel function \Hq~\x)\ 2 ~ -. That the phase 9 
is constant implies that a coherence behavior appears in 
the system, i. e. the localized state is a phase-coherent 
state, as previously discussed 34| . It is a unique feature of 
wave localization, and has also been shown to be related 
to electronic localization (e. g. Ref. H^). 

Spatial behavior of localized waves Following |36|, a 
general consideration of the spatial behavior of localized 
waves is possible. Consider a wave transmitted in a ran- 
dom medium. The transport equation for the total en- 
ergy intensity /, i. e. < \E\ 2 >, may be intuitively written 
as 



dl 

— = -al, 
dx 



(11) 



where a represents decay along the path traversed. After 
penetrating into the random medium, the wave will be 
scattered by random inhomogeneities. As a result, the 
wave coherence starts to decrease, yielding the way to 
incoherence. Extinction of the coherent intensity Iq, i- e. 
< E > | 2 , is described by 



die 

-r- = -lie, 
dx 



(12) 



with the attenuation constant 7. Eqs. (|11|) and (|12f> lead 
to the exponential solutions 



I(x) = I(Q)e~ 



and I c (x) = l{0)e~ 



(13) 



When there is no absorption, the decay constant a 
is expected to vanish and the total intensity will then be 
constant along the propagation path. Then, the coherent 
energy gradually decreases due to random scattering and 
transforms to the diffusive energy, while the sum of the 
two forms of energy remains a constant. This scenario, 
however, changes when localization occurs. Even without 
absorption, the total intensity can be localized near the 
interface due to multiple scattering. When this happens, 
a does not vanish. The transport of the total intensity 
may be still described by Eq. Ullfl . and the inverse of a 
would then refer to the localization length. 

The above perceptual description may be illustrated 
by Fig. |3J Without or with little absorption and with 
no localization, the energy propagation is anticipated to 
follow the behavior depicted in (A). When the localiza- 
tion comes in function, the wave will be trapped within 
an e-folding distance from the penetration, as prescribed 
by (B). In the non-localization case, the diffusive in- 
tensity increases steadily as more and more scattering 
occurs, complying with the Milne diffusion. In the lo- 
calized state, the diffusion energy increases initially. It 
will be eventually stopped by the interference of multiple 
scattering waves. Issues may be raised with respect to 
whether this apprehended image is supported by actual 
situations. In the rest, we will inspect this problem. 
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In deriving these equations, the boundary condition was 
used; it states that 1(0) — lc(0) as no scattering has 
been incurred yet at the interface. According to energy 
conservation, the incoherent intensity Ijj (diffusive) is 
thus 

I D (x) =I(x)-I c (x). (14) 



FIG. 3: Picturised non-localization and localization effects. 

Comparing Fig. [3] (A) and (B), we know that one of 
the key differences between non-localization and localiza- 
tion features lies in the difference in the behavior of the 
diffusive waves. However, the eventual task of inferring 
localization is to find whether the total energy decays 
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with the sample size. 

C. Numerical results 

Unless otherwise noted, the following parameters are 
used in the numerical simulation: the non-dimensional 
damping rate, &o/ w o = 0.001 and the interaction cou- 
pling, b = 0.001. The filling factor (/?) varies from 2.25 
to 25; the filling factor is defined as the number of dipolcs 
per unit area. But without notification, the filling factor 
is taken as 6.25. The number of random configurations 
for averaging is taken in such a way that the convergency 
is assured. In the calculation, we scale all lengths by a 
length D such that koD — 1, and frequency by ujq. In 
this way, the frequency always enters as k/k. We find 
that all the results shown below are only dependent on 
parameters b, bo/u>o, and the ratio u/too or equivalcntly 
fc/fcrj. Such a simple scaling property may facilitate de- 
signing experiments. In the numerical computation, wc 
take c = 1 for convenience. The total wave at a spatial 
point is scaled as T(r) = E(r)/Eo(r), with Eq being the 
direct wave from the source, so that the trivial geometric 
spreading effect is naturally removed. 

First, we plot the frequency response of the transmis- 
sion in the scenario (B) of Fig. The results are shown in 
Fig. 01 Here we see that there is a narrow window within 
which the transmission is highly inhibited, implying a 
strong localization effect. It is also clear that when the 
sample size is increased, the inhibition increases. Com- 
paring the results for the two filling factors, we know that 
the strong inhibition regime increases with filling factor. 
In the next computations, we will focus on frequencies 
within the strong inhibition region. 



Here we see clearly shown that for the three frequencies 
within the strong inhibition region, the energy is spatially 
confined near the transmitting source, and, as expected, 
the energy seems to decrease nearly exponentially along 
any radial direction. Meanwhile, the system reveals an 
in-phase phenomenon: nearly all the phase vectors of 
the dipoles point to the same direction, exactly opposite 
to the phase vector of the source which is denoted by 
the black arrow. The picture represented by Fig. [5] fully 
complies with the general description of the coherence in 
localization stated above. 
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FIG. 5: The phase diagram for the two dimensional phase 
vectors defined in the text; the phase of the source is assumed 
to be zero. Each vector is located at the site of the dipole; 
thus the locations of the phase vectors also denote the random 
distribution of the dipoles. Bottom: The spatial distribution 
of energy (~ \T\ 2 ). Three frequencies are chosen, the sample 
size is 8 x 8. 
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FIG. 4: Transmission versus frequency for two filling factors 
and two sample sizes. 

Now we consider the phase and the spatial distribu- 
tion of energy of the system. To describe the phase 
behavior of the system, we assign a unit phase vector, 
Ui = cos6ie x + sin9ie y to the oscillation phase 6i of the 
dipoles. Here e x and e v are unit vectors in the x and 
y directions respectively. These phase vectors are repre- 
sented by a phase diagram in the x — y plane with the 
phase vector tt, being located at the dipole to which the 
phase 9i is associated. The results are depicted for three 
frequencies by Fig. [S] 



We also note from Fig. [5] that near the sample bound- 
ary, the phase vectors start to point to different direc- 
tions. This is because the numerical simulation is carried 
out for a finite sample size. For a finite system, the energy 
can leak out at the boundary, resulting in disorientation 
of the phase vectors. When enlarging the sample size 
by adding more dipoles while keeping the averaged dis- 
tance between dipoles fixed, the area showing the phase 
coherence will increase accordingly. 
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FIG. 6: Comparison of the total transmission at three fre- 
quencies for the two scenarios shown by Fig. [5] The x-labels 
'Distance' and 'Size' refer to scenario (A) and (B) respectively. 
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The results of Fig. [5] are encouraging, as they are a 
strong indication of localization. In the following, we will 
further explore the features of localization. In Fig. El we 
compare the transmission results for the two scenarios 
from Fig. [5] Here is shown that although there is a slight 
difference in the transmitted strength, overall speaking 
the spatial decay features are nearly identical, signified by 
the match of the decaying slopes indicated in the figure. 
Though suspected to be true previously, such a match 
is important, and to the best of our knowledge this is 
the first that has been ever shown for EM waves. It 
supports directly that the scenario (B) can also be used 
to infer localization effects, facilitating measurements of 
localization; in the original conjecture, it was scenario 
(A) that has been suggested for discerning localization. 
In the rest of simulation, we will adopt scenario (B) in 
Fig.0 

The bottom panel of Fig. El indicates that the level of 
spatial localization of energy varies for the three frequen- 
cies. To quantify the localization in Fig. El we plot the 
total energy as a function of the sample size. The results 
are presented by Fig.0 Here, the numerical data are fit- 
ted with the least squares method and the fitted curves 
are shown by the solid lines; the unnoticeable deviation 
from the lines reflects the fluctuation due to the random 
distribution. Two ways of averaging are adopted. One 
is the traditional way in which the logged total trans- 
mission is averaged, while the other is to take log of the 
averaged total energy. These can be referred to from the 
y-axis labels. It shows that after removing the spreading 
factor, the data can be fitted by e~ r ^. From the slope 
of the solid lines, the localization lengths £, which are 
the inverse of the slopes can be estimated. Here is shown 
that the decaying slopes from the two averaging methods 
are very close, an encouraging fact. It is also indicated by 
the decreasing slopes with frequency that the localization 
effect decreases as the frequency increases. 



consideration of localization. Due to scattering and lo- 
calization, the coherent waves decrease with the sample 
size. The diffusive wave increases initially as more and 
more scattering occurs, then reaches a peak and starts 
to decay due to the localization effect, in agreement with 
the description of Fig.|3{B). The results show that in the 
present system, the diffusive portion in the total energy is 
much smaller than the coherent portion, indicating that 
the mean free path is very small. When plotted in the 
log scale, we have found that the total and coherent en- 
ergies decays exponentially with the distance, while the 
diffusive energy increases initially, then starts to decay 
exponentially. It is worth noting that for the three fre- 
quencies considered, the total energy does not reveal any 
behavior similar to the diffusive waves, in contrast to the 
previous theoretical conjecture^- The theory predicts 
that the total energy would follow the behavior of dif- 
fusive waves until when the sample size is larger than 
the localization length. One explanation is that in the 
present system, the scattering is too significant so that 
the diffusive portion never dominates. A search for the 
possible match between the simulation and theory in cer- 
tain conditions is still undergoing. 
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FIG. 7: The transmission versus the sample size for three fre- 
quencies. The slopes can be used to estimate the localization 
length. 

With the fixed filling factor of 6.25, we have also in- 
vestigated the spatial variations of the total, coherent 
and the diffusive energies for three frequencies discussed. 
The results are presented in Fig. El Here we see that 
the results are in accordance with the previous general 



FIG. 8: Behaviors of the total, coherent, and the diffusive 
energies as a function of sample size in the scenario described 
byFig.HB). 




(B) 






-A- wi<o = 1 .003 




co/ro = 1 .005 




to/(i> = 1 .007 



1.002 1.004 1.006 1.008 1.01 

Frequency(o>/oi ) 



10 15 20 25 
Filling Factor)!!) 



FIG. 9: Localization length as a function of (A) frequency 
and (B) filling factor. The sample size is 8 x 8 for (A). 

In Fig. El we plot the localization length versus fre- 
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quency and filling factor separately. It is shown that 
with the fixed filling factor, referring to Fig. EfA), the 
localization length increases with frequency within the 
frequency regime considered. With a fixed frequency, the 
localization length tends to decreases, meaning increas- 
ing localization effects, as the filling factor increases. 

Fig. I1UI shows the transmission versus frequency for 
various coupling constants and damping rates. From this 
figure, we observe the following. (1) The increasing cou- 
pling strength leads to a wider strong inhibition region, 
but shallower localization valley. (2) When increasing the 
coupling strength, a prominent resonance peak appears 
below the natural frequency ujq. And the peak moves 
toward lower frequencies as the strength increases, a fea- 
ture also appears in the acoustic system [2^. (3) Overall 
speaking, the increasing damping rate degrades the lo- 
calization level, and tends to abolish the resonance peak. 
Also it seems to widen the strong localization region at 
the lower frequency side. 



the localization depends on a number of parameters in- 
cluding frequency, filling factor, and damping rate. The 
spatial behavior of the total, coherent and diffusive waves 
is also explored, and found to comply with a physical in- 
tuitive picture. A phase diagram characterizing localiza- 
tion is presented, in a gree ment with previous investiga- 
tions on other systems [34] . 
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IV. SUMMARY 

In this article, the localization features in a simple elec- 
tromagnetic system are investigated in detail. Some gen- 
eral properties of the localization phenomenon are elab- 
orated. For certain ranges of frequencies, strongly local- 
ized electromagnetic waves have been observed in such a 
simple but realistic disordered system. It is shown that 



FIG. 10: Transmission versus frequency for various cou- 
pling strengths and damping factors: (A) various coupling 
strengths with damping rate 0.001; (B) various damping rates 
with coupling constant 0.001. The sample size is 4 x 4. 
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